Dispersion Coefficients by a Field-Theoretic Renormalization of Fluid Mechanics 
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We consider subtle correlations in the scattering of fluid by randomly placed obstacles, which have 
been suggested to lead to a diverging dispersion coefficient at long times for high Peclet numbers, in 
contrast to finite mean-field predictions. We develop a new master equation description of the fluid 
mechanics that incorporates the physically relevant fluctuations, and we treat those fluctuations by 
a renormalization group procedure. We find a finite dispersion coefficient at low volume fraction of 
disorder and high Peclet numbers. 
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PACS numbers: 47.10+g,05.60.Cd,47.55.Mh 

Dispersion of fluid particles by flow through random 
media is important in problems ranging from pollutant 
transport through soil to enhanced oil recovery to chem- 
ical reactor design. Recently, it has been suggested that 
the dispersion coefficient may diverge logarithmically at 
long times [[!]. Such a result, which is not inconsistent 
with experimental data on flow through packed beds, 
would have significant implications. This prediction, 
however, is in disagreement with mean-field results that 
predict a finite dispersion coefficient Accounting for 
fluctuations with a renormalization group approach, we 
find the dispersion coefficient should be finite. 

There are four sources of randomness not accounted for 
in the mean-field treatment of fluid dispersion. The nu- 
merical results suggest a logarithmic divergence in three 
dimensions, i.e. that the upper critical dimension may be 
three, and the mean-field theory, which lacks the random- 
ness, would not be expected to detect such a divergence. 
The first source of randomness is the random motion of 
the fluid particles. Since the Navier-Stokes equations are 
the starting point for the mean-field theory, all micro- 
scopic fluctuations of the fluid particles have been sup- 
pressed. A second source of randomness is the random lo- 
cations of the obstacles to the fluid flow. The mean-field 
theory is essentially a renormalized single-obstacle the- 
ory. As such, all correlations around an obstacle beyond 
the Brinkman screening length are ignored. The third 
source of randomness present in the numerical simula- 
tions, but not present in the mean-field theory, are meso- 
scopic fluctuations naturally present in the lattice Boltz- 
mann simulation method. Imperfections in the pseudo- 
random number generator are a fourth and final source 
of randomness present in the simulations. Using a field- 
theoretic renormalization group approach, we take into 
account the first three physical sources of randomness. 
We find that these sources of randomness do not modify 
the mean-field prediction: the dispersion coefficient re- 
mains finite at long times and high Peclet numbers. The 
Peclet number is given by Pe = v x R/D{, where v x is the 
average velocity of the fluid, R is the average radius of 



the blocking obstacles, and Df is the diffusion coefficient 
of the pure fluid. 

Our goal is to write a field theoretic representation of 
the Navier-Stokes equation: 

dtVt + ILfc Ej VjdjVk = v\l 2 v l + fa (1) 

where v = /i/p is the kinematic viscosity, and = 
n^-Ffc — dkP)/p is the total body force on the fluid. 
The presence of the projection operator 11^ (k) = 5^ — 
kikk/k 2 in these formulas ensures that the incompress- 
ibility condition V • v = is maintained || . The Fourier 
transform is defined by /(k) = J dx/(x) exp(ik ■ x). 

We first write a master equation model of fluid me- 
chanics. Generalizing previous treatments from one di- 
mension Q] to d — 3 dimensions, we consider the fluid 
on a lattice, with Na "velocity" particles on lattice site 
A. Since the velocity is a vector, we denote the number 
of particles contributing to the i th component of the ve- 
locity as N x . The state of the system is defined by the 
configuration of velocities in the system, and the master 
equation relates how the probability of any given config- 
uration, P({iV A }), changes in time 

d t P = AP (2) 

where A is a linear operator that specifies the dynam- 
ics of the fluid. Defining the identity, /, shift operators, 
E l x P{N l x ) = P(N{ + 1), E X "P(N X ) = P(N{ - 1), and 
a number operator, N{P(N{) = N{P(N{), the dynam- 
ics of Eq. H are given by diffusion, convection, and force 
terms, A = A d + A c + A f : 
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Here h is the lattice spacing, and 8u is the contribution 
of one velocity particle to the velocity. The sum is over 
lattice sites for A and over the three dimensions for jkl. 
The vector gives the lattice site displaced by one unit 
in the k th direction. Expanding this master equation for 
small lattice spacing, one finds that it reproduces the 
Navier-Stokes equation in this limit. 

We map this master equation onto a field theory using 
the coherent states representation ^| . Negative veloci- 
ties are treated as anti-particles. In this mapping, we find 
the correspondences 6*(r) <-> E\ , bi(r) «-> E\N x /h d , 
and b*(r)bi(r) <-> N x /h d . As is typical, we set b*(r) = 
b~i(r) + 1. We find the following action: 



S = 



dt bi(-k, t)[d t + vk 2 + S(t)}b l (k,t) 
-iXi f f dt (2^) d (5(k 1 +k 2 +k 3 ) 

x fc 2 .Sfe(k 1 ,t)S^(k 2 ,t)S J ± (k 3 ,i) 

I I dt {2n) d 5(k 1 + k 2 + k 3 + k 4 ) 

x fc 2 . t k (k x , t)bi (k 2 , t)b^ (k 3 , t)bf (k 4 , t) 
dtfA(0,t) (4) 



Here the notation J k stands for J dk/(27r) d , the integrals 
over time are from t = to some large time t = tf, 
and the summation convention is implied. This action 
is written in terms of the divergence-free part of the ve- 
locity, 6^~(k) = n i fc(k)6fc(k). For convenience in later 
calculations, we have included a curl-free component in 
the quadratic terms (Feynman gauge). Initially, A< = 1. 

This field theory for fluid mechanics differs from the 
traditional one J?], |) in that random fluctuations of the 
fluid are incorporated by the presence of the A 2 term. 
Using this action, standard results from fluid mechanics 
are reproduced. For example, a calculation of the long- 
time tails in the velocity-velocity correlation function in 
two-dimensions is in accord with the known result [01: 
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(5) 

The scaling of the Kolmogorov energy cascade and the 
Richardson separation law are also reproduced for the 
common statistical model of turbulence 
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E{k) ~ (const) £T 5/3 
r 2 (t) ~ (const)* 3 



(6) 



In both of these calculations, the A 2 term does not con- 
tribute. This is in contrast to recent work on reaction- 
diffusion systems, where fluctuations contribute at low 
density |5|, |6| |ll| . Since the fluid density is always finite 



in the present case, the density fluctuations captured by 
A 2 prove to be irrelevant. 

With these preliminaries taken care of, we now turn 
to the problem of flow through porous media. The main 
effect of the porous media is to exclude the fluid from 
certain fixed obstacles. The velocity is zero within and 
on the surface of these obstacles. Within the context 
of the master equation, if an obstacle is at site A, the 
velocity there must be zero, N x = 0. This fixes the 
lattice spacing as h w R, where R is the characteristic 
size of the blocking particles. A relatively singular model 
for the obstacles sets the average fluid velocity at the 
obstacle sites to zero. With this model, the obstacles 
generate the following additional term in the action: 
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dt b\b\ 
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here P[n^] = (1 — <j))5 n 



where cf> is the prob- 



ability of a site being occupied by an obstacle, i.e. the 
volume fraction of obstacles. This condition is difficult 
to implement directly. The condition [J dt b\b\] 2 < a 2 , 
where a 2 \/(AnvR) 2 is a constant, effectively imple- 
ments this condition in the long-time limit and is much 
more convenient to implement. That is, we use a Gaus- 
sian to regularize the delta function 
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Expanding to the lowest relevant order in the fields, we 
find the additional term in the action to be 



5" = (const) 



d d x / dtdt' 



2a 2 h d 

x 6 J (x,t)o l (x,t)6 4 (x,i')o l (x,i') (9) 



We set 71 = cf>/(2a 2 h d ) for later convenience. We note 
that this form looks like a contribution arising from ran- 
dom, imaginary, velocity-dependent point forces acting 
on the fluid. 

A key aspect of the flow through porous media is the 
net fluid flow. To accommodate this, we shift the fields by 
the average values. We consider the flow to be along the 
x direction, and so set bi = 5i^ x v x + b[ and then rewrite 
the action in terms of b\. Suppressing the primes, we find 



S = 



'k 



dt bi(-k,t)[d t + vk 2 - ik x v x + 5{t)]k(k,t) 
dt (27r) d ,5(ki+k2+k 3 ) 
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Xk h b k (k!,f)^(k 2 ,t)&4-(k3,i) 

f f dtdt 1 (27r) d (k! + k 2 + k 3 + k 4 ) 

«/kik2k3k4 J 
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x6, (ki,t)^(k a ,f)6, (k 3 ,t')bl(k 4 ,t') 
dtdt' (2 7 r) d (k 1 +k 2 +k 3 ) 



kik 2 k 3 



+73 



xt> x (ki,t)S£(k 2 ,t)S I (k 3 ,0 

y y dfatt' b x {-Kt)b x (k,t') - j & fSt(o,t) 

(10) 

where 72 = 2v x ji and 73 = V x ji. We have eliminated 
the A2 terms, as they do not contribute in the renormal- 
ization procedure. To incorporate the random obstacles, 
we have used the replica trick but have suppressed 
these details that do not enter in a one-loop calculation. 

We now apply the renormalization group procedure 
to this action. From power counting, the upper criti- 
cal dimension for this theory is d c — 3. The frictional 
force on the fluid due to the obstacles will generate a 
mass term, however, that eventually renders the theory 
finite in any dimension. We use the momentum shell 
procedure, where fields on a shell of differential width 
dl are integrated out, Ae~ dl < k < A, where A a ir/li 
is the cutoff, and I is the flow parameter. As usual, we 
rescale time by the dynamical exponent t' = te~ zdl , dis- 
tance perpendicular to the flow direction by k' ± = k±e dl , 
and distance along the flow direction by the dilation ex- 
ponent k' x = k x e vdl We make use of the average 

(C(ki)^(k2)> = (27r) d 5(k 1 +k 2 )n ij (k 2 )G (k 2 ) where 
Go(k) = \/(vk 2 — ik x v x ). We perform a one-loop cal- 
culation, valid to first order in disorder strength and all 
orders in the parameter Ai. A large number of integrals, 
over one-hundred, appear in this renormalization proce- 
dure, and details will be presented elsewhere. A typical 
integral would be one such as 
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on shell J -00 l/ ±^l_ + V x^ X 



shell 



u±v x kj 



4irv x v± 



dl 



(11) 



We have made use of the fact that the integral over k x 
converges, and so it can be taken off the shell jyj), and 
we have used the fact that v x flows to zero and can be 
ignored in these calculations, i.e. we replace vk 2 by v\k\ 
in the propagator, where k 2 = k 2 + k\. 

When computing the renormalization flows, some ad- 
ditional terms are generated. We present the de- 
tails only for the rc-fields, as these are most signifi- 
cant. A term such as 7iu / klk2k3k4 / dtdt' (27r) d (ki + 



k 2 +k 3 +k 4 )b x {k u t)b^(k 2 ,t)b x (k 3 ,t')b^(k 4 ,t) is gen- 
erated. In addition, a mass term is generated, 
rrixj^jdt b x (—k,t)b x (k,t). Defining the dimension- 
less couplings g x = 71 / '(4irv±v x ) , g xu = 7iu/(47r2/_i_u r ), 
92 = 72/(871x1 and g 3 = 7 3 /(47w_ L v 3 ,), we find the 
following flow equations: 

din g 3 



dl 
dm 92 
dl 

d In g x 
dl 

d In g xu 
dl 

d In 
dl 

dm T 



(2z-d+l 



2g X u 



dl 

din v x 

dl 
d In v\ 



= (2z - d + 1 - rj) - 2g x 

+Xi x {g x +g xn )g 3 /g 2 
= (2z - d + 1 - rj) - 2g x + 2X lx g 2 - X 2 lx g 3 

= (2z-d+l-r))- 2g xu ~ Ag x + 4X lx g 2 
+2Xi x g 2 g x /g xu - >4 x 939x/g X u - 2X\ x g 3 
= (z - 77) - 2g x - 2g xu + 4X lx g 2 - 2Xl x g 3 

= zm x + 2g x vj_A 2 + 2g xu v±A 2 - 2X lx g 2 v±A 2 

din/ 



= (z ~ V) - 4Xix92 + 3A lx 33, 
(U = ~ 2) - 2X lx g 2 - ^X\ x g 3 



dl 



(12) 



At long times, we find g x (l) ~ g 2 (l) ~ g 3 (l) ~ l/(g 1 + 
21). We find the unphysical term g xn {l) should vanish. 
We also find X\ x ~ A?/(l + 2lgo). From the flow equa- 
tion for the viscosity, we find the dynamical exponent 
is z = 2 + 0(l~ 2 ). From the flow equation for the ve- 
locity, we find the dilation exponent is r\ — 2 + 0{l~ 2 ). 
The mass quickly reaches the asymptotic form m x (l) ~ 
m*e 21 whatever the non-universal, A-dependent terms in 
dm x /dl. In the absence of the mass term, e.g. for imag- 
inary, velocity-dependent forces with zero average value, 
these flow equations generate interesting scaling behav- 
ior. In the present case, the generated mass stops the 
flow at a characteristic value of the flow parameter, I*. 

We determine the value of the generated mass by a 
force balance argument. From the requirement that 
(b x ) = 0, we find m x (l)v x = f(l), and conclude that 
fo = m*v x . From fluid mechanics, it is known that 
the force density is given by Stokes' Law at low vol- 
ume fraction of obstacles and moderate flow rates ||: 
fo = 6ttRi/±v x <Po/(4:TtR 3 /3), where R is the radius of 
the spherical obstacles. Equating these two results, 
we find m*/v±_ = 94> /(2R 2 ). We further calculate 



the correlation length in the fluid £ x 



.21 



Ml)/m x (l)} 1/2 



e 2l [v e~ 



2 7(m*e 2i )] 1/2 : 
1/2 



(13) 



which is exactly the Brinkman screening length ||. At 
high flow rates, empirical expressions for the drag force 
are available that can be used to identify m* |L4| . 
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We now address the issue of the dispersion coefficient. 
The dispersion coefficient D x = L dt'[{v x (x, t)v x (x + 
icv x t' , t + t')) — v 2 } is a measure of the "dispersive" trans- 
port induced by the fluid. We can calculate this quan- 
tity directly from the field theory by perturbation theory. 
We find D X {1) = g 3 (l)v 2 /m(l). Using the approximate 
value of a 2 and the determined form of m(l), we find 
D x (l) = e 4nRv x /9. We match the flow equations to 
this perturbation theory when the flow equations begin 
to loose their validity, m x {l*) « v±A 2 , which gives the 
characteristic flow parameter I* = i ln(Vj_A 2 /m*). We 
calculate the dispersion coefficient from matching |ll|] as 
D* — e 21 D x (l*). We find the renormalized dispersion 
coefficient to be 
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D* 



Air 
~9 



Rv x 



(14) 



This dispersion coefficient is due solely to mechanical dis- 
persion of the fluid, and the form is in exact agreement 
with the mean-field result of Koch and Brady Q . The 
dispersion coefficient is proportional to the Peclet num- 
ber and independent of the volume fraction in the limit of 
high Peclet number and low volume fraction. Only the 
numerical prefactor, here given approximately, is non- 
universal. Interestingly, the asymptotic value of the dis- 
persion coefficient is independent of I*, as long as the 
mass has been generated and l*g « 1, which is the case 
for small volume fractions or large Peclet numbers. The 
independence of the observable on degree of renormal- 
ization is common when mean-field theory is essentially 
exact, as is the case here. 

In the dilute limit, the dispersion coefficient grows as 
shown in figure |l]. The dispersion coefficient reaches its 
asymptotic value at a time on the order of to. If time 
is nondimensionalized as r = tv x /R, this characteristic 
time is given by r = Pe 1 / 3 0, which is the time scale 
for diffusive transport across the boundary layer. 

The suggestion that dispersion coefficients diverge log- 
arithmically in three dimensions at large Peclet numbers 
and finite volume fractions of disorder H, then, is not 
borne out by the present calculations. Such a divergence 
would show up as a I* dependence of the predicted disper- 
sion coefficient, which is absent in our calculations. The 
proposal of a diverging dispersion coefficient was based 
upon simulation data, and these simulation data may 
not have sampled the full diffusive boundary layer trans- 
port that occurs on a time scale that grows as Pe 1 / 3 . 
While transport at high volume fractions is, in principle, 
not accessible by our renormalization group calculations, 
physical aspects of the expected behavior are present in 
our result. In particular, curing of the logarithmic diver- 
gence in the theory by growth of a mass term due to the 
frictional drag force exerted by the disorder particles on 
the fluid seems generic. This frictional drag force would 
seem to prevent divergences at any volume fraction of 
disorder. Indeed, the higher the volume fraction of disor- 
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FIG. 1: Shown is the behavior of the dispersion coefficient 
at low volume fraction of disorder. The dispersion coefficient 
reaches its asymptotic value at a dimensionless time of tq = 
Pe 1 ' 3 . 



der, the greater the expected accuracy of the mean-field 
result. Our result is in agreement with the mean-field 
theory as soon as the frictional mass has been generated. 

The renormalization group treatment of dispersion 
for a particle-based model shows that mean-field theory 
misses no essential physics at low volume fractions. 
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